Non-linear simulations of combustion instabilities with a 
quasi- ID Navier-Stokes code 



Nils Erland L. Haugen, 0yvind Lang0rgen and Sigurd Sannan 

SINTEF Energy Research, N 0-7 465 Trondheim, Norway 



Abstract 

As lean premixed combustion systems are more susceptible to combustion 
instabilities than non-premixed systems, there is an increasing demand for 
improved numerical design tools that can predict the occurrence of combus- 
tion instabilities with high accuracy. The inherent non-linearities in com- 
bustion instabilities can be of crucial importance, and we here propose an 
approach in which the one- dimensional Navier-Stokes and scalar transport 
equations are solved for geometries of variable cross-section. The focus is 
on attached flames, and for this purpose a new phenomenological model for 
the unsteady heat release from a flame front is introduced. In the attached 
flame method (AFM) the heat release occurs over the full length of the flame. 
The non-linear code with the use of the AFM approach is validated against 
results from an experimental study of thermoacoustic instabilities in oxy-fuel 
flames by Ditaranto and Hals [Combustion and Flame, 146, 493-512 (2006)]. 
The numerical simulations are in accordance with the experimental measure- 
ments and both the frequencies and the amplitudes of the resonant acoustic 
pressure modes are reproduced with good accuracy. 

Keywords: Combustion instabilities, thermoacoustics, numerics 
1. Introduction 

New and stricter emission regulations, particularly on nitrogen oxides, 
have led to the adoption of lean premixed combustion as the primary tech- 
nology for low-emission power generation from gaseous fuels in stationary 
gas turbines. Homogeneous premixing of the fuel and the oxidizer, combined 
with an ultra-lean operation mode, provide for lower combustion temper- 
atures and a drastic reduction in NO x formation, accompanied by lower 
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emission of soot and CO. Experience has shown, however, that lean pre- 
mixed combustion systems are more susceptible to combustion instabilities 
than non-premixed systems. The coupling between unsteady heat release 
and acoustic pressure oscillations can lead to self-excited oscillations that 
cause unacceptable levels of noise and moreover tend to reduce combustion 
efficiency These thermoacoustic oscillations also have a detrimental effect 
on the combustor equipment that limits the component lifetime and, in the 
worst case, can result in system failure due to structural damage. 

The thermoacoustic instabilities involve a feedback cycle in which the 
heat release rate and acoustic pressure fluctuations are coupled. A stability 
criterion was deduced by Lord Rayleigh jlj, which states that the pressure 
wave will be amplified and develop an instability if the heat release fluctu- 
ations are in phase with the acoustic pressure oscillations. Oppositely, the 
system is stable if the heat release rate and acoustic pressure fluctuations 
are out of phase. In general, there are various types of combustion instabil- 
ities due to the different physical mechanisms that can drive the instability. 
The underlying driving mechanism depends on the combustor geometry, the 
burner and flame types, and the set-up of the fuel and oxidizer feed lines. 
Thus, an instability due to variation in the equivalence ratio may arise if the 
acoustics is able to interact with the upstream flow all the way to the indi- 
vidual feed lines. Another type of instability occurs if the acoustic pressure 
oscillations affect the incoming velocity of the fuel-oxidizer mixture into the 
flame front. Reviews of the acoustically coupled combustion instabilities, as 
well as other instability mechanisms that do not couple directly to the system 
acoustics, are given by Candel |i[ and McManus et al. j^j. 

Combustion instabilities and the interactions between several different 
physical phenomena such as variations in the flow velocities, acoustic pres- 
sure fluctuations, and heat release are inherently complex by nature. Tra- 
ditionally, combustion system designers have relied on experiments in order 
to obtain vital information for the design of combustors where instabilities 
can be avoided or controlled. However, performing a series of large-scale 
experiments, in some cases under high pressure to create gas turbine operat- 
ing conditions, is very costly. For this reason the development of numerical 
design tools based on Computational Fluid Dynamics (CFD) has become an 
important endeavor in recent years in order to complement experiments in 
the design process. The complexity of the instability mechanisms and their 
interactions, however, places severe demands on both the physical modeling 
and the numerical analysis of current CFD techniques. Hence, the devel- 



2 



opment of CFD as a reliable design and analysis tool with sufficiently high 
prediction accuracy with regard to combustion instabilities is an extremely 
challenging task. 

The most accurate method to study turbulent combusting flows numeri- 
cally is by Direct Numerical Simulation (DNS), where the flow field is solved 
directly from the Navier-Stokes equations pj. Due to the huge computational 
cost, however, the use of DNS is restricted to simplified problems and turbu- 
lent flows with relatively low Reynolds number. For high Reynolds number 
flows in complex geometries, as often encountered in industrial applications, 
DNS is therefore prohibitively expensive and is likely to be beyond reach for 
quite some time. A computationally less demanding approach is Large Eddy 
Simulation (LES), in which the large turbulent eddies of the flow are com- 
puted explicitly, while the smaller eddies are not resolved but modeled using 
a subgrid scale model. LES is well suited to the description of unsteady dy- 
namical phenomena, and there has been a growing interest in LES in recent 
years |5| . But even LES is limited due to expensive computational costs, and 
thus 3D LES calculations are non-practical for use on an everyday basis in 
industrial applications. 

As a consequence of the complexity of instability mechanisms, much of the 
modeling work on combustion instabilities has focused on simplified models 
to make problems tractable. The development of linear models, in which the 
Navier-Stokes equations and scalar transport equations are linearized, has 
been a very useful approach capable of predicting combustion instabilities at 
a qualitative level [3j. Thus linear acoustic models can predict, at least to 
some extent, the frequencies of resonant modes and their growth rate during 
a phase of exponential growth. A current practice in the modeling of unsta- 
ble combustion systems is to apply a network model where the geometry of 
the combustor is modeled by a network of acoustic elements and a simplified 
form of the pressure equation is solved. The acoustic elements of a network 
model, also called multiports, correspond to various components of the sys- 
tem, e.g., the air or fuel supply, the transition between two regions of different 
cross-section, the outlet of the combustor, or the flame itself 0, Q, 0, 0, 10 . 
Mathematically, each multiport in the system is represented by a transfer 
function. 

Although the predictive scope of linear models is limited, the use of lin- 
ear models for active control of combustion instabilities has been widespread 
[3], However, linear models are not able to predict the amplitudes of the 
resonant modes of the statistically stationary state, and hence the dominat- 
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ing instabilities of the system cannot be distinguished from the less significant 
ones. Also, the couplings between the instability mechanisms may not be ac- 
curately represented in linear models and there has therefore been a growing 



interest in active control based on non-linear models [121 ] . When the non- 
linear Navier-Stokes and scalar transport equations are solved, it is possible 
to predict the shape and the level of the acoustic frequency spectrum both 
during the exponential growth phase and for the statistically stationary state. 
This should make for a more accurate and efficient active control of the com- 
bustion instabilities. In a non-linear real space description it is furthermore 
possible to apply more reliable models for the heat release, as discussed in 
Sec. 14.21 An alternative approach to solving the non-linear real space equa- 
tions is to solve a non-linear wave equation in the frequency domain Lj, |l4 . 

In this paper we have chosen an approach where the non-linear Navier- 
Stokes and scalar transport equations are solved for a combusting flow in real 
space and time, as in a DNS or an LES, but in a reduced one-dimensional 
description. Previous work on real-space ID non-linear CFD simulations of 



combustion instabilities include the work of Polifke et al. [l(| and Rook 15 



Polifke et al. [10[ considered thermoacoustic oscillations in a straight duct 
by introducing a one- dimensional heat release model for a heat source placed 
in the duct. From time-dependent simulations of the dynamical behavior 
of the heat source, the authors were able to obtain the transfer matrix of 
the heat source which again could be investigated in a linear network model. 



Rook [15| applied one- dimensional CFD to study the acoustical behavior 



of a burner-stabilized flame using a pressure correction method to simu- 
late the flame on a ceramic foam burner. We here follow the approach of 



variable cross-section area introduced by Cohen et al. [16L and also used by 



Prasad and Feng [18|, [19| , where the one-dimensional equations are written for 
variable-area geometries in order to simulate three-dimensional geometries. 
The reduction of the governing equations to ID is valid for relatively simple 
geometries for which the flow can be assumed to be quasi-one-dimensional, 
and where the wavelengths are sufficiently large so that the wave motion is 
well approximated by a plane wave. In this framework we propose a new 
phenomenological model for the unsteady heat release from a flame front. 
The flame model, here termed the attached flame method (AFM), gives a 
ID realization of the flame front. Using AFM it is demonstrated that the ID 
simulations capture both the exponential growth and the nonlinear statisti- 
cally stationary state of the acoustic oscillations at a computational cost far 
below that of a corresponding DNS or an LES. 
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In Section [2] we present the governing equations, and make them one- 
dimensional by introducing a variable cross-section. Boundary conditions is 
discussed in Section [31 while the new flame model is described in Section HI 
Finally, the code is verified in Section |5] by comparing the simulation results 
with experimental results from a study of combustion instabilities in oxy-fuel 



flames by Ditaranto and Hals [20 



2. The governing equations 

The governing equations for a turbulent combusting flow are based on the 
conservation of mass, momentum and energy, and the transport equations 
for species mass fractions. The conservation of mass is represented by the 
continuity equation 

f + V-pu = 0, (1) 
while the conservation of momentum gives the Navier-Stokes equations 

du 

p— + pu- Vu = -Vp + V-t (2) 

where p is the density, u is the velocity vector, p is the pressure, and r is 
the viscous stress tensor. For a Newtonian fluid the stress tensor is given 
by r = 2pS, where p is the dynamic viscosity, and S is the traceless strain 
tensor with components 

1 / dm duj _ 2 du k \ 
ij 2 {<),■, ' <),-, :C J Or,)' ' ; ' 

when we use the Einstein summing convention. The conservation equation 
for the energy can be rewritten in terms of the temperature equation 

dT 1 
p— + pu- VT= - -pV-u + V- (AVT) + g, + g c ), (4) 
at c v 

where T is the temperature, c v is the specific heat at constant volume, A is 
thermal conductivity, q v = 2/iS 2 is the viscous heating source term, and q c is 



The bulk viscosity has been set to zero by Stokes' hypothesis. 
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the heat release rate from combustion. The equations for the mass fractions 
Y k of the species k can be written as 



QYi 

p^- + p u • VF fc = V • ( P DVY k ) +pu k ; k = l,...,N s , (5) 

where D is the mass diffusivity, u k is the chemical reaction rate of species k, 
and N s is the number of species. The above set of equations are closed via 
the ideal gas equation of state 

p = prT, (6) 

where r is the gas constant of the mixture. The mixture gas constant is given 
by r = R/m, where R = 8.31 kJ/ (kmol K) is the universal gas constant, and 
m is the mean molar mass. The speed of sound in the gas mixture is 



c 



yfyrT, (7) 



where 7 = c p /c v is the specific heat ratio, with c p the specific heat constant 
pressure. 

2.1. The variable cross-section ID approximation 

In a variable cross-section geometry the combusting flow can be solved by 
a quasi-one- dimensional treatment. The reduced governing equations in one 
dimension are presented here. In a quasi-lD description the traceless strain 
tensor is reduced to 

/2 1 1\ du 

s = diag (r-3--3k' < 8) 

From this it follows that the viscous force is 

v_v,^s>^(g + g^), m 

where p is the dynamic viscosity. Furthermore, the viscous heating reduces 
to 

=!"(£)"■ (10) 

It should be noted that the factor 4/3 in the above equations is due to the 
ID approximation. 
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The quasi- ID continuity equation for a variable cross-section geometry 
can be written as 

dp = 1 dpAu 

dt A dx ' 1 ' 

where A is the cross-sectional area. It follows that the one-dimensional equa- 
tions for the momentum, the temperature and the species mass fractions are 
given by 

du du 1 dp 4 / d 2 u ldudp\ Ff )W 
dt dx p dx 3 y dx 2 pdxdxj p 

dT dT 1 / p duA d , #Z\ \ 

"aT = ~ u IT + — 7~5 — + 7r\ X Jr) +Qv + q c + qv,w (13) 

dt dx pc v \ A dx dx K dx J J 

^ 9Yk dY k Id ( n dY k , 

u^— + —^-[pD-—] +u kl (14) 



dt dx pdx^ dx 

where Ff )W represents the viscous force from the walls, and q V)W is the corre- 
sponding viscous heating. These terms are added to the system since the vis- 
cous force in Eq. fl9]) contributes to damping in the streamwise direction only, 
i.e., wall effects do not naturally appear in the one dimensional equations. 



For more information on the viscous force Ff tW , see Appendix |Appendix A 
If the viscous wall effects are negligible, Ff )W and q v , w in Eqs. f fT2|) and f fT3|) 
can be set to zero. 

Numerically, the set of equations (JTTJ-flHJ) are solved using an explicit 
solver. The spatial discretization is sixth-order finite difference, while third- 
order Runge-Kutta is used for the time stepping. 

3. Boundary conditions 

3.1. Open boundaries 

We assume that the time varying pressure and velocity field in the system 
can be decomposed as 

P = Po+P 

u = Uq + u' (15) 

where subscript denote the mean part while primes denote the fluctuating 
part. In addition, the density fluctuations p' are given by 

P = Po + p\ (16) 
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where p is the average density. The acoustic pressure fluctuations are de- 
scribed by plane harmonic waves given by 



p' = p e - wt = (A + e ikx + A~e~ ikx ) e~ iut , (17) 

where k = 27r/A is the wave number, u is the angular frequency, and A + 
and A~ denote the amplitudes of the right-moving and left-moving pressure 
waves, respectively. Similarly, the velocity fluctuations are given by 



u > = far** 



— ( A + e lkx - A~e~ ikx ) e~ iwt . (18) 



PqCq 

In a long duct with an open exit at x = 0, the impedance at the exit is 

z = ! = poCo (i+^f) = PoCo (iTl) ' (19) 

where the reflection coefficient R is defined by R = —A~/A + . Rearranging 
the above equation we get 

R = E2E°-Z. (20) 
PoC + Z 

We notice that if Z = there is perfect reflection and R = 1. If Z = poCo, 
corresponding to the characteristic impedance of the medium, R = and 
there is no reflection at the exit. Since p'/u' = p/u = Z, we may use the 
equation of state (jSJ), together with Eq. f|T9l) . to write 

1 - (u-u). (21) 
Solving this equation with respect to the temperature, we obtain 

r + -10 (£!)). (22) 

Thus, at an open boundary the temperature is determined by the expression 
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3.1.1. Reflection coefficients at open boundaries 

We consider acoustic oscillations in long ducts. If the duct is a circular 
pipe, the characteristic Helmholtz number is defined as H n = ka, where a 
is the radius of the pipe. Thus, the Helmholtz number is small if the pipe 
radius is small compared to the acoustic wavelength. It is known that for 
low Helmholtz numbers the dominating acoustic modes within a pipe will be 
reflected from an open end, i.e., most of the acoustic energy will not leave 
the pipe jilj . This is normally also the condition for combustion instabilities 
to appear. By assuming that H n <C 1 and applying conservation of mass, 
it has been shown that the absolute value of the reflection coefficient at the 



open end of an unflanged pipe is [22 

\R\ = 1 - ^(ka) 2 ; ka < 0.2. (23) 

For large acoustic velocities u, non-linear effects at the exit will no longer 
be negligible. The dominating non-linearity being the vortex shedding at the 
sharp bends of the exit. Peters et al. [22) find the acoustic power absorbed 
by vortex shedding, and non-dimensionalized by P nor m = ^pu 3 7ra 2 , to be 

Cortex = ^ = /3SrV 3 (24) 



norm 



for the acoustic Strouhal number Sr ac = ^ 3> 1, and 

-^Vortex = TTT (25) 

for Sr ac C 1, In the above q = 2 for an unflanged pipe and we have set 
(3 = 0.5 in the simulations. 

If Re(Z) is the real part of the impedance Z, the power lost at the outlet 
due to radiation from the pipe exit can be found by using Eq. ( fl9|) and 
Eq. (J23D to be 

P rad = 7ra 2 / = na 2 ^Re(Z)\u\ 2 , (26) 
where I is the intensity. This gives 

Pr * ad = = \ kaSTac - (27) 
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Combining the above equations gives the following expression for the reflec- 
tion coefficient function of the power loss 



R=f -p, (28) 

CO I p* ' v 1 

u ' loss 

when P{ oss = P* OTtcx + -P r * ad is the sum of the acoustic losses through radiation 
and vortex shedding. 

3.2. Closed boundaries 

Closed boundaries are understood as acoustically closed boundaries that 
reflect acoustic waves. Thus, closed boundaries can either be closed or open 
with respect to the mass flow. If there is no inflow (mass flow) across the 
acoustically closed boundary the velocity u, together with the derivatives |^ 
and are set to zero. On the other hand, if there is inflow across the 
acoustically closed boundary, the mass flow is given a constant value (the 
velocity u has a non-zero constant value). 

An acoustically closed boundary with inflow can be thought of as wall 
with several small holes, e.g., a porous plate. The holes are so small that 
they do not affect the reflecting abilities of the wall. However, they are large 
enough so that the mass flow entering the domain is significant and is allowed 
to cross the boundary. 

By setting ^ = the closed boundary is adiabatic. This is not necessarily 
precisely correct, but the error introduced by this assumption is expected to 
be of minor importance. 



4. The flame model 

In this section we obtain an expression for the combustion heat release 
term q c in Eq. ( IT3|) . We begin by describing the well-known n-r formulation 
and show how it can be integrated into a Navier-Stokes solver. We then 
introduce a new phenomenological flame model which we refer to as the 
attached flame method (AFM). 

4-1. The n-r model for unsteady heat release 

As a first approximation we apply the n-r model, which provides a global 
description of the unsteady heat release rate associated with combustion 
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instabilities. The heat release rate (energy per unit time) in the n-r model 
is given by 

Q' = ^-u(x,t-r), (29) 
7-1 

where A is the cross-section of the duct, n is the interaction index determining 
the coupling between the velocity and heat release fluctuations, and r is the 
time lag between these fluctuations. The heat release rate per volume is 

€ = £.= ^u( X , t -r), (30) 

where Lf is the flame length, and we have made use of the relation ([7j). 

The n-r model was originally designed for linear wave-equation models 
for which the main focus has been on the heat release fluctuations and their 
coupling to the flow-field perturbations. For an implementation of the n-r 
formulation in a CFD code based on the Navier-Stokes equations, the mean 
part of the heat release needs to be included as well. Thus, the total heat 
release rate per volume (the mean part plus the fluctuation) can be written 

(ivy T \ 
h c + —j—u(x,t - r) J , (31) 

where h c is a constant heat source. Both the constant h c and the unsteady 
heat release are in this approach non-zero only within the flame. 

4-2. The attached flame method (AFM) 

The n-r formulation is primarily designed to be used with linear wave- 
equation models in the frequency domain. Although generalizations to non- 
linear applications have been made, the n — r model is not optimal for use in 
a non-linear Navier-Stokes solver where the governing equations are solved 
in real space and time. Thus, in the n-r model the heat release is limited to 
a point source. In addition, the time lag r must be accounted for. 

More detailed analyses of unsteady heat release from flame fronts have 
been based on studies of the dynamics and shape of anchored premixed flames 



23l l24j . Fleifil et al. [25[ used a kinematic model to calculate the transfer 
function in the linear regime between heat release and upstream velocity 
oscillation of a premixed flame stabilized on the rim of a tube. The kinematic 



approach was later extended by Schuller et al. |26| by including convective 



effects of the flow modulations propagating upstream of the flame. In this 
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paper we present a new phenomenological model for the unsteady heat release 
from an attached flame. The model, denoted the attached flame method 



(AFM), is similar to the analytical flame front model |25|, l26( in the sense 
that a real-time differential equation is solved to capture the kinematics of the 
thermoacoustic instabilities. However, while the kinematic model is based 
on the G-equation approach to determine the location of the flame front, 
in AFM the flame position is obtained directly from the equations for the 
species mass fractions. The basic idea is here to project the flame, which is 
essentially two-dimensional, into a one- dimensional description. 

One advantage of the AFM approach, compared to the n — r formulation, 
is that for laminar flames no free parameters such as n and r need to be 
determined. Another advantage is that the heat release occurs over the full 
length of the real flame, and not just from a single point. In the AFM 
formulation, the full length Lf of the flame is defined as the distance from 
the inlet, where the flame is anchored, to the point where all the fuel is 
burned. The flame length Lf is therefore a dynamical variable that changes 
with time as the flow conditions change. 

The laminar flame speed is known to vary slightly with the position in a 
flame. Thus, the laminar flame speed will in general be different at the base 
of the flame, in the flame tip, and in the main body of the flame. The spatial 
variations in the flame speed are believed to have only minor effects on the 
flame front, however, and have been neglected in the following. A known 
constant laminar flame speed has been assumed in the current formulation, 
although this can easily be changed if the effects of a variable flame speed 
are expected to be significant. 

For a combusting flow in a long duct we assume that the flame front can 
be represented by piecewise straight lines, as shown in Figure [TJ The total 
flame surface AAi within a ID grid cell is then given by 

AAi = 2AH i( l, (32) 

where d is the depth of the ID grid cell, i.e., the size of the grid cell in 
the direction perpendicular to the paper plane. The flame front is here the 
interface between the fresh and the burned gas within a given ID grid cell. 
Since the approach is one-dimensional only, the flame front is not resolved 
by the ID code but it is evident that the distance from the centerline to the 
flame front, denoted hi, is proportional to the mean mass fraction Yf ue i of the 
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Figure 1: In this duct of height 2h there is a flame holder at the position Xf in the 
streamwise direction. The flame front is represented by piecewise straight lines with total 
length H = A-ffj, and the length of the flame is Lf. Every ID grid cell has length Ax. 



fuel in the ID grid cell. This can be written as 

Yfael 



hi 
h 



(33) 



fuel, inlet 



where ifu e i,iniet is the mass fraction of the fuel at the inlet. We note that the 
limiting cases of hi = h upstream of the flame anchor at Xf and hi = down- 
stream of the flame tip are recovered from the above equation. Differentiating 
Eq. (I3"3"|) we obtain 

dhi dY iue \ h 



dx 



dx Yf 



fuel, inlet 



which, when setting Ahi 



Ax^, gives 



Ax 2 + Ahj = Ax 



\ 



1 + 



fuel 



Y 



fuel, inlet 



dx 



(35) 



In the above derivation we have used the fact that the mass fraction of the 
fuel outside the flame cone is zero (burned gas), while the mass fraction of 
the fuel inside the flame cone equals the mass fraction of the fuel at the inlet 
(fresh gas). 
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Since the flame front is always moving into the fresh gas, the reaction 
rate of the fuel within the flame is given by 



„ AAiY {uel 

,inlet*J Z/JT /nn\ 

R= ^ , (36) 

where Sl is the laminar flame speed, and AV is the volume of the grid cell. 
The reaction rate of the fuel is then given by 

{0 x < Xf 

- 25Lyf t in ' ct/ V 1 + (^^) 2 ^<*<^ + Lt (37) 
x > Xf + Lf. 

where fx > 1 is a constant accounting for the fact that the real flame speed 
might be larger than the laminar flame speed due to turbulence. In fact, 
fx, quantifying the turbulence in the flame, is the only closure needed in the 
model. In the current work the focus is on laminar problems, however, and 
the model is here validated for laminar or weakly turbulent cases only, for 
which fx — 1. 

By construction the AFM does not allow for more than two flame fronts 
for any given downstream position. That is; a given flame segment can not 
have an angle of more than 90 degrees to the wall. This means that the 
model will not be able to describe e.g. vortex roll up, and is therefore most 
suited to compact flames. This is a natural consequence of the model being 
one dimensional. 

For a chemical reacting system with Ns species, the general equation for 
a chemical reaction can be written 

N S N S 

E -» E (38) 

k=l k=l 

where A k symbolizes the chemical species k, and v' k and v' k ' are the stoichio- 
metric coefficients of species k on the reactant and product side, respectively. 
The chemical reaction rate u k for such a system may be expressed as 

w fc = K' - ^D-ftfuei-^, (39) 

"ifuel 

where m k and mf ue i is the molar mass of species k and of the fuel, respectively. 
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The energy release within a grid cell is determined by the reaction rate 
.Rfuei and the lower heating value h L of the unburned mixture. Thus, the 
heat release term q c in Eq. ffl3l) is defined by 

q c = R lucl h L p. (40) 

The expressions (139]) and (|40|) are used to close the set of governing equa- 
tions dHD-(H3D. 

4-3. The secondary grid 

In many applications of premixed combustion the geometry is such that 
a narrow slot leads into a wider combustion chamber. This is visualized in 
Figure [2] where a slot of height h s \ ot leads into a combustion chamber of height 
/i com b. In such a geometry it is no longer correct to use the same convective 
velocity for the species convection in the jet as in the mean flow. The rea- 
son is that the mean convective velocity of the flow within the combustion 
chamber is lower than the convective velocity of the fresh fuel/air jet entering 
from the slot. Hence, we here solve the velocity evolution of the jet using 
a separate secondary subgrid. The jet entering the combustion chamber is 
visualized by the thick dashed lines in Figure [2j which then also represent 
the upper and lower boundaries of the secondary grid within the chamber. 
Note, however, that since the simulation tool is formulated in ID only, no 
boundary conditions are required at the dashed lines. In the streamwise di- 
rection the secondary grid is bounded by the coordinates x sec l and x sec>2 , as 
shown in Figure [2j In this subgrid domain the governing equations (lli p - ()14p 
are solved for a constant cross-section A and with no reactions, i.e., Uk and 
q c are both zero. This gives the convective velocity u conv to be used in the 
convective term of Eq. (|12[) . The secondary grid is thus applied to only a 
small fraction of the entire computational domain covered by the primary 
grid, which means that special care must be taken concerning boundary con- 
ditions. At the secondary subgrid inlet x SCCj i the boundary values are chosen 
as the corresponding instantaneous values of the primary grid at the same lo- 
cation. For the subgrid outlet at x sec 2 the non-reflecting boundary condition 
given by R = in Eq. (1201) is used. 

It should be pointed out that the matching of the secondary and the 
primary grid solutions in the above manner does not have any significant 
impact on the final results as long as both the inlet and the outlet of the 
secondary grid are outside the flame and a non-reflecting boundary condition 
is used at the secondary grid outlet. 
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Figure 2: The thick solid lines represent the geometry of a slot leading into a combustion 
chamber. The thick dashed line corresponds to the upper and lower boundary of the 
secondary grid within the chamber. The coordinates x SCCi i and x SCC! 2 define the boundaries 
of the secondary grid in the streamwise direction. 



5. Validation of the quasi-lD code 

The non-linear quasi-lD code is here validated against analytical, ex- 
perimental and numerical results. We first compare with analytical results 
obtained in a simplified setup of a flame in a duct. The non-linear code is 
then validated against results from an oxy-fuel study in a lab-scale test rig, 
both for inert and reactive flows. 

5.1. Straight duct 

We here consider a simplified case of a flame in a straight duct for which 
an analytical solution is known when the flame is described by the n-r model. 
For this comparison case the AFM model is not applied since a direct com- 
parison between the AFM and the n-r model is not straightforward, given 
that various choices of n and r may lead to very different physical solutions. 
Thus, the quasi- ID numerical results are here compared with the analytical 
solution given by Poinsot and Veynante [fj, with the n-r model applied in 
both approaches. 
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The simulated duct is assumed to be a straight pipe of length 2a = 1 m, 
with the flame located in the center of the pipe. For the analytical treatment 
the heat source is a point source, while for the non-linear solver the flame 
has a length of 8 cm in the pipe direction. The spatial extension of the flame 
is required since any gradient or object in a spatial code must be resolved 
by at least a few grid points. The set-up is such that air enters at room 
temperature at the pipe inlet. The constant heat source h c given in Eq. f l3Tj) 
is chosen to be zero such that the mean temperature at the outlet equals the 
inlet temperature. The interaction index n between the velocity and the heat 
release is set to n = 0.25. We are here primarily interested in the imaginary 
part of the frequency, which for this simplified case is given by 

Jm (^i) = ~ sm(^r), (41) 
where c is the speed of sound, and 

2tic 

Uj = 2wf j = (l + 2j) — . (42) 

od 

In this linear approach the amplitude of a resonant frequency Uj is given 
by Aj = Bj exp(— icujt), such that the amplitude Aj grows exponentially 
when the imaginary part of Uj is positive. If the imaginary part of Uj is 
negative, on the other hand, the amplitude of this particular frequency will 
decay exponentially. In Table I5TT1 the imaginary parts of the first four u/s 
are shown for three different values of r. For r = 6 x 10~ 4 s it is evident 
that oo i = 1639 Hz and U3 = 3824 Hz are the unstable frequency modes. 
In the left plot of Figure [3] the positions of these two modes are marked by 
large arrows, while the smaller arrows correspond to the decaying resonant 
modes. The solid line of the plot shows the energy spectrum obtained from 
the numerical quasi-lD simulation, and it is seen that the excited modes 
agree well with the unstable modes as given by the analytical solution. In 
the last column of Table [57T1 the simulated growth rate is also found to be in 
good agreement with the theoretical values. For r = 15 x 10~ 4 s we note that 
the unstable frequency modes are given by u\ = 1639 Hz and u 2 = 2731 Hz. 
In the middle plot of Figure [3] we observe that the amplitude of U\ is much 
larger than the amplitude of u 2 , despite the fact that the imaginary part of 
u 2 is larger than that of U\. This is due to that the initial perturbations 
are such that Bi » B 2 , and that the amplitude A 2 , while still in the linear 
regime of the simulation, did not have enough time to catch up with A\. 
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This is also reflected by the fact that the simulated growth rate corresponds 
to the imaginary part of Finally, for r = 27 x 10~ 4 s we find that there 
are no unstable frequency modes and that all the amplitudes are decaying. 
The resonant modes are nevertheless recovered in the right plot of Figure |3l 
but the amplitudes are very weak (note the scale on the ordinate axis). For 
this case it should be noted that due to quite large uncertainties in the 
determination of the decay rate, the growth rate presented in Table 15.11 has 
relatively large error bars. 




10 3 10 3 10 3 

co [Hz] co [Hz] co [Hz] 



Figure 3: Energy spectra for simulations with r — 6 x 10 4 s (left), t = 15 x 10 4 s 
(middle) r = 27 x 10~ 4 s (right). 



r[10" 4 s] 


Im(u ) 


Im(ui) 


Im(u 2 ) 


Im(u 3 ) 


Simulated growth rate 


6 


-2.2 


5.8 


-6.9 


5.2 


5.2 


15 


-5.1 


4.4 


5.7 


-3.6 


4.3 


27 


-6.9 


-6.6 


-6.2 


-5.4 


-5.6 



Table 1: Growth rate of unstable frequencies for different r's. The angular frequencies are 
w =546 Hz, wi=1639 Hz, w 2 =2731 Hz, and w 3 =3824 Hz. 



5.2. Sudden expansion burner 

In this section we validate the non-linear quasi- ID cod e ag ainst experi- 



mental results from the oxy-fuel study of Ditaranto and Hals [20|]. A schematic 



view of the geometry of the lab-scale combustion rig is shown in Figure HI 
The set-up consists of a 80 cm long premixing section, followed by a 4 cm 
long flame arrestor and a 10 cm long slot leading into the 47.9 cm long com- 
bustion chamber. The premixing and combustor sections are square channels 
with cross-section s x s, where s = 5.4 cm is the inner width of the sections. 
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For the flame arrestor and the slot the cross-sections are s x d arr . and sxd 
respectively, where d arr , = 2 cm and d s i ot = 0.5 cm. 



I — gas inlets 



slot ) 




N 2 purge | 
port 3 « I > Q port 4 




mixing plate flow homogenizer 



sudden expansion 



Figure 4: Schematic view of the geometry of the oxy-fuel combustion rig used in the 
experimental study of Ditaranto and Hals [2fjJ . 



5.2.1. The cold rig 

In order to validate the non-reactive part of the non-linear code, acoustic 
simulations were performed for the geometry of the oxy-fuel rig described in 



20i |. with no flame and only air present in the rig. Experimentally, in the cold 
rig case a loudspeaker was placed at the upstream end of the premixer, i.e., 
at the inlet at the left-end side of the set-up shown in Figure HI The acoustic 
frequency response of the loudspeaker was measured with four microphones, 
and for this procedure there was no flow in the system. The loudspeaker 
additionally produced white noise for which the spectral distribution is not 
known. The numerical simulations can therefore not reproduce the ampli- 
tudes of the resonant modes, but merely aim at reproducing the resonant 
frequencies. 

The velocity spectrum corresponding to the statistically steady state so- 
lution of the non-linear code is plotted in Figure Also shown in this figure 
are the experimentally measured values of the resonant frequencies, along 
with the resonant frequencies obtained with the use of the linear code. We 
note that there is generally a good match between the measured resonant 
modes and those obtained by the numerical simulations. The resonant peak 
at ~ 25 Hz corresponds to the low frequency 1/4 mode of the entire system. 
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Figure 5: The cold flow velocity spectrum from microphone port 3 of the rig illustrated 
in Figure |4] The solid curve shows the spectrum obtained from the non-linear code. The 
upper arrows indicate the measured resonant frequencies, while the lower (smaller) arrows 
indicate the resonances obtained from the linear code. 



5.2.2. Oxy-fuel combustion 

We here validate the full quasi- ID non-linear code using the AFM formu- 
lation against an experimental study of combustion instabilities in oxy-fuel 
flames in a sudden expansion test rig. In the oxy-fuel study of Ditaranto 



and Hals [20| four different combustion instability regimes, referred to as 
regions, were distinguished. In Figure |6] are shown typical spectral distribu- 
tions obtained in the oxy-fuel experiment for two of these regions. The shown 
pressure spectra were recorded at microphone port 3 illustrated in Figure HJ 
In what is referred to as Region 3 the flame front is attached to the slot of 
the rig at all times. For Region 2 the oxy-fuel flame study exhibited two 
instability patterns; one for which the flame is not attached to the slot but 
follows the formation of periodic large vortices, and another corresponding 
to a hysteresis phenomenon for which there is a combination of an attached 
flame branch and vortex shedding. 

We observe from the spectral distributions of Figure Elthat the instability 
frequencies are quite different for the shown cases. Thus, for the attached 
flame of Region 3 the 1/4 mode from the combustion chamber is at a signifi- 
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cantly lower frequency (~ 230 Hz) than the corresponding modes (~ 300 Hz) 
for the cases in Region 2. According to Ditaranto and Hals [2oj, the instability 



Region 3 (Re=2000) 

Region 2 (Re=2000) 




Frequency [Hz] 



Figure 6: The spectral distributions of pressure oscillations in the oxy-fuel rig as measured 
by Ditaranto and Hals 20] . The pressure spectra were recorded at microphone port 3 of 
the rig in two different instability regions. In Region 3 (thick solid line) the flame is 
attached to the slot. In Region 2 there are two instability patterns; one (thin solid line) 
displaying periodic vortex shedding and another (thin dashed line) being an intermediate 
case with both an attached flame branch and the formation of periodic vortices. The left 
arrow in the figure points to the low- frequency 1/4 mode of the entire system, the three 
next arrows point to the frequency peaks of the combustor 1/4 mode, the two upper arrows 
point to the peaks of the premixer 1/2 mode, while the two arrows to the right point to 
the frequency peaks of the 1/1 mode of the premixer. 

frequency for the attached flame case is in fact a combination of the premixer 
1/2 mode and the combustion chamber 1/4 mode. The premixer 1/2 modes 
for the cases of Region 2 both have frequencies at ~ 200 Hz, although the 
mode corresponding to the hysteresis case is not dominant and difficult to 
observe in Figure |6j Similarly, the premixer 1/1 modes for the cases of Re- 
gion 2 have frequencies at ~ 400 Hz, while the corresponding mode of the 
attached flame case is at ~ 460 Hz. The main difference between the two 
instability cases of Region 2 is that the premixer 1/2 mode is most amplified 
when the flame follows the shedded vortices entirely, while the combustor 
1/4 mode is dominant when the flame is partially attached to the slot. 



21 



The combustion chamber temperature was not measured by Ditaranto 



and Hals [20| and we therefore use the frequency peak of the ground mode of 
the combustor to estimate the average temperature. From Figure |6] we note 
that the 1/4 mode of the combustion chamber is at a frequency fo ~ 230 Hz 
for the attached flame. We first compute the speed of sound and recall that 
the combustor has length 47.9 cm. For a duct with an open end we also know 
that the acoustics depends on an end correction to the duct length. For an 
open end unhanged pipe, Davies et al. [27| have obtained an empirical fit to 
the end correction given by 

l/a = 0.6133 -0.1168 (ka) 2 ; ka < 0.5, (43) 

where a is the radius of the pipe. For the square combustion duct we set 
a = s/y/rr = 3.0 cm, corresponding to a pipe of equal cross-section as the 
duct. This gives a length correction of I ~ 1.8 cm for the ground mode. 
Adding this to the combustor length, we find the acoustic length of the 
combustor to be / ac0 ustic = (0.479 + 0.018) m = 0.497 m. With f = c /A and 
A = 4/ acoustic for the 1/4 wave mode, this gives the average speed of sound in 
the combustion chamber 

Co = 4/o^acoustk = 458 m/s. (44) 

Using Eq. ([7]), the average temperature in the combustion chamber then 
becomes 

c 2 

T = -2- = 720 K, (45) 
7r 

where 7 = 1.24 and r = 235 J/(kg K) for the given gas mixture. The average 
temperature is thus much lower than the adiabatic flame temperature at 
about 2400 K (when assuming a combustor inlet temperature of 400 K). 
In the following we therefore allow for sufficient cooling in the combustion 
chamber such that the mean temperature downstream of the flame is 720 K. 

As discussed in Section 14.21 the AFM model is not able to describe vor- 
tex roll up, but is designed to describe attached compact flames. For the 
validation of the quasi- ID code a numerical simulation of the case for which 
the flame was attached to the slot has therefore been performed. The case 
is typical of the instability regime of Region 3, and defined by a Reynolds 
number Re = 2000, an equivalence ratio $ = 0.9, and a volume fraction of 
42% 2 in the 2 /C0 2 oxidant. The corresponding pressure spectrum is 
shown in Figure O In Ditaranto and Hals 20| the spectral distribution of 
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the pressure oscillations is shown in Fig. 7c. From Fig. 1 in 20] the laminar 
flame speed is Si = 0.47 m/s in the case of 42% 2 in the oxidant. 

In the experimental study the flow rate of fuel and oxidant in the at- 
tached flame case was 1.3 x 10~ 4 kg/s and 1.65 x 10~ 3 kg/s, respectively, 
corresponding to a mean flow velocity in the premixer of 0.42 m/s. With 
Re = 2000 the attached flame is in the laminar flame regime. Hence, the 
parameter f T in Eq. (J3"7|) is set to = 1, i-e., the flame speed is equal to 
the laminar flame speed. 

The observed thermoacoustic instabilities in the flame experiments are 
governed by statistically-stationary limit cycles in which acoustic pressure 
variations lead to fluctuations in the flow velocity and heat release. The heat 
release fluctuations, on the other hand, feed back to the acoustic modes at 
the same frequency. For the oxy-fuel study the limit cycles are controlled 
by a saturation in the heat release caused by acoustic losses at the outlet. 
This can be observed in the left graph of Figure [71 where the envelope of the 
pressure amplitudes at microphone port 3 of the rig is shown as a function of 
time for various simulations. We note that when there are no acoustic losses, 
i.e., when the reflection coefficient R — 1, the instability grows to infinity. 
Taking linear acoustic effects into account, we obtain from Eqs. ( 123]) and ( 14*4]) 
a reflection coefficient of R = 0.995 at the dominating frequency fo = 230 Hz. 
The corresponding instability growth is shown by the thick dashed curve in 
Figure [7] and displays a slightly smaller growth rate than for the R = 1 case. 
By including both non-linear and radiative losses, the reflection coefficient 
can be calculated dynamically from Eq. f ]28]) . The resulting instability growth 
is shown by the solid curve denoted "Non- linear" in Figure [71 If, in addition, 
viscous damping from the walls of the flame arrestor is taken into account, 
i.e., the term Ff w in Eq. ( Il2p is non-zero, the growth rate is given by the 
fairly similar dashed curve denoted "Damp.+non-lin.". (For more details 
on Ff >w , see Appendix Appendix A ). Finally, the dashed-dotted curve in 
Figure [7] shows the pressure envelope for a simulation for which R = 0.96. 

From the right graph of Figure [7] we observe that the thick dashed curve 
and the solid (thin) curve have very similar slopes for small pressure varia- 
tions, while for larger pressure amplitudes the non-linear losses become more 
and more important and the corresponding solid curve levels out at a much 
smaller amplitude than the dashed curve. Comparing the simulation results 
when viscous damping effects were taken into account (dashed curve) with 
no damping effects implemented (solid curve), we observe that the viscous 
damping only has little effect on the pressure envelope except for a smaller 
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growth rate initially. The simulation using a reflection coefficient of R = 0.96 
was done as a comparison case and produced a slightly larger pressure am- 
plitude than those indicated for the solid and dashed curves. But the initial 
growth rate produced by the simpler model produced an initial growth rate 
that was smaller. From these observations we conclude that it is of crucial 
importance to the numerical simulations that non-linear damping through 
vortex shedding at the combustor outlet is included. The viscous damping 
at the walls turns out to be of less importance for this specific set-up. One 
issue to be kept in mind, however, is that the equations (121)1 and ([2"7]) for 
non-linear and radiative losses, respectively, were deduced for an open end 
circular pipe. For the application of a square duct considered here, additional 
losses due e.g. to the duct corners might therefore have an impact. 

In Figure [8] the pressure spectra obtained experimentally (thick solid 
curve) and by numerical simulations (thin dashed and solid curves) are 
shown. We observe that all the main resonant frequencies of the experi- 
ment are well matched by the simulations. In addition, there is a fairly good 
agreement between the experimental and numerical values of the levels of 
the resonant peaks of the pressure spectra. However, the levels of the curves 
in between the resonant modes are lower for the numerical simulations than 
what was obtained in the experiment. This may be explained by the fact 
that experimental peaks are normally broader than their numerical counter- 
parts, but the possibility that the difference is due to the one-dimensional 
approximation can clearly also not be excluded. It is also interesting to note 
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that the peaks of the two numerical pressure spectra are well aligned, except 
for a small difference for the higher frequency modes. 




Figure 8: Experimental and numerical pressure spectra in the combustion chamber at 
microphone port 3. 



6. Conclusion 

In this paper we first present the quasi- ID Navier-Stokes and scalar trans- 
port equations for a variable cross-section duct. The temperature and species 
mass equations are closed by expressions for the heat release and the reac- 
tion rates obtained from the AFM formulation. The AFM approach has 
been introduced as a new phenomenological flame model in real space and 
time for describing attached flames in a one-dimensional geometry of variable 
cross-section. With the use of the secondary subgrid it is shown that two- 
dimensional features of a jet entering the combustion chamber is accounted 
for. 

The quasi- ID code is non-linear and gives the time evolution of real space 
quantities. This is in contrast to conventional linear wave-equation models 
that only give the growth rate of the unstable resonant frequencies, and under 
the assumption that a linear representation is adequate. Thus, the quasi-lD 
code is capable of reproducing the non-linear saturation of the acoustic oscil- 
lations, the so-called limit cycle. In order to achieve the correct magnitude 
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of the limit-cycle oscillations, it is crucial to account for the acoustic losses 
at the open end(s) of a duct. This is done by including the non-linear effects 
due to vortex shedding at the sharp bends at the duct exit, in addition to 
the losses due to acoustic radiation from the open end exit. 

The non-linear code has been validated by first comparing results from 
using the simplified n-r heat release model with results obtained with the n-r 
model in a linear wave-equation solver. The code was then validated against 
the oxy-fuel study of Ditaranto and Hals [20| in a sudden expansion burner. 
With the application of the AFM heat release formulation the numerical sim- 
ulations were able to reproduce the resonant frequencies of the the acoustic 
pressure spectrum of an oscillating attached flame with high accuracy. In 
addition, the levels of the resonant peaks were reproduced quite well. From 
these findings we conclude that the quasi-lD non-linear Navier-Stokes solver 
with the AFM formulation is a promising tool for further studies of combus- 
tion instabilities in a variety of cases, including jets in variable cross-section 
ducts or in co-flows. 
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Appendix A. Viscous damping 

The acoustic Strouhal number is defined by Sr ac = ua/u, where u is the 
angular frequency of the acoustic oscillations, a is the radius of the pipe, 
and u is the amplitude of the acoustic velocity oscillations. If the acoustic 
Strouhal number is very small, the boundary layer develops in a time much 
smaller than the acoustic period. In that case it is a good approximation to 



assume that the boundary layer is always developed. Following White 28 
the viscous force from the walls Ff tW then is 



f,u> 



-r w SpAx, (A.l) 
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where t w is the wall shear stress, Sp is the duct perimeter, and Ax is the 



length of a grid cell. The wall shear stress can be expressed as [28 



(A.2) 



where fjj is the Darcy friction factor. The Darcy friction factor can be 
approximated with the expression 

64/prrf,! 



D 



Re Dh 



(A.3) 



for laminar flows. For turbulent flows fjj can be obtained from the relation 

fprof,2^D h 



2.0 log 



0.8. 



(A.4) 



VJd 

The Reynolds number is based on the hydraulic diameter = 4A/Sp, 
where A is the cross- section of the duct, such that Re^ = uD^/p. In the 
above equations, / pro f,i and / pro f,2 are factors whose values are determined by 
the shape of the duct. Thus, for a flow in a circular pipe / pro f,i = / pr of,2 = 1- 
For a flow between two parallel plates we have / pro f,i = 3/2 and /p ro f,2 = 0.64. 
It should be noted that the expression (1A.2I) is valid only for fully devel- 



oped flows, see for instance the work by Allam and Abom 29] for further 
details. This limitation is most prominent for pipes with large cross-sections, 
i.e., when the Strouhal number Sr ac is not much smaller than 1. In a flame 
trap where the viscous damping in general is the largest, however, the given 
expressions should be relatively good due to the very small cross-sections of 
the holes. For a more applicable expression for acoustic losses in pipe flows, 
see the paper by Disselhorst and Wijngaarden 
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Following Peters et al. [22 



for small Helmholtz numbers ka and large 
shear numbers the damping coefficient due to the viscous boundary layer is 
given by 



a = — 

Co 



V2Sh 



1 + 



v/PT J Sh 2 



1 



2Pr 



(A.5) 



where Sh = ayu/v is the shear number and Pr is the Prandtl number. The 
above expression for a has been established under the assumption of no mean 
flow in the system. This gives a viscous damping defined as 

AP' ~ exp(-ax), (A.6) 

where AP' is the decrease in the acoustic pressure fluctuations, and x is the 
distance traveled by the acoustic waves. 
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